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Abstract 

Three high order shock-capturing schemes are compared for large eddy simulations (LES) 
of temporally evolving mixing layers (TML) for different convective Mach numbers (M c ) 
ranging from the quasi-incompressible regime to highly compressible supersonic regime. 
The considered high order schemes are fifth-order WENO (WEN05), seventh-order WENO 
(WEN07) and the associated eighth-order central spatial base scheme with the dissipative 
portion of WEN07 as a nonlinear post-processing filter step (WEN07fi). This high order 
nonlinear filter method (H.C. Yee and B. Sjogreen, Proceedings of ICOSAHOM09, June 
22-26, 2009, Trondheim, Norway) is designed for accurate and efficient simulations of 
shock-free compressible turbulence, turbulence with shocklets and turbulence with strong 
shocks with minimum tuning of scheme parameters. The LES results by WEN07fi using 
the same scheme parameter agree well with experimental results of Barone et al. (2006), 
and published direct numerical simulations (DNS) work of Rogers & Moser (1994) and 
Pantano & Sarkar (2002), whereas results by WEN05 and WEN07 compare poorly with 
experimental data and DNS computations. 

Key words: High order numerical methods, Numerical methods for turbulence with 
shocks, DNS, LES, Mixing layer 


1. Introduction 

Part of the inaccuracy in direct numerical simulations (DNS) and large eddy simulations 
(LES) of turbulent flow using standard high order shock-capturing schemes is due to the fact 
that this type of computation involves long time integrations. Standard stability and accu- 
racy theories in numerical analysis are not applicable to long time wave propagations and/or 
long time integrations [38]. The original construction of modern shock-capturing schemes 
was developed for rapidly developing unsteady shock interactions and short time integra- 
tions. Any numerical dissipation inherent in the scheme, even for high resolution shock- 
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capturing schemes that maintain their high order accuracy in smooth regions (e.g., fifth- 
or seventh-order WENO schemes (WEN05 and WEN07)), will be compounded over long 
time integration leading to smearing of turbulence fluctuations to un-recognizable forms. 

In compressible turbulent combustion/nonequilibrium flows, the constructions of nu- 
merical schemes for (a) stable and accurate simulation of turbulence with strong shocks, and 
(b) obtaining correct propagation speed of discontinuities for stiff reacting terms on “coarse 
grids” share one important ingredient - minimization of numerical dissipation while main- 
taining numerical stability. Here “coarse grids” means standard mesh density requirement 
for accurate simulation of typical non-reacting flows. This dual requirement to achieve 
both numerical stability and accuracy with zero or minimal use of numerical dissipation is 
most often conflicting for existing schemes that were designed for non-reacting flows. In 
addition to the minimization of numerical dissipation while maintaining numerical stabil- 
ity in compressible turbulence with strong shock, Yee & Sjogreen, Yee and Yee & Sweby 
[46, 47, 43, 42] discussed a general framework for the design of such schemes. Yee & 
Sjogreen [51], Sjogreen & Yee [37, 36] and Wei et al. [39, 40] and references cited therein 
present their recent progress on the subject. In [55], a short overview of this recent progress 
is given. The discussion addresses three separate yet interwoven types of numerical chal- 
lenges for high speed turbulent reacting flows containing discontinuities. This paper is 
confined to the comparison of three specific high order shock-capturing methods on tur- 
bulent mixing for non-reacting flows. The study for turbulent mixing for reacting flows is 
planned. 

2. Recent Progress in Numerical Methods for Turbulence with Strong Shocks 

The current trends in the containment of numerical dissipation in DNS and LES of tur- 
bulence with shocks are summarized in Yee & Sjogreen and Yee et al. [51, 50, 54]. See 
the cited references for details on these current trends. Before presenting the improved 
filter schemes and their application to the temporally evolving mixing layers (TML), the 
key ingredients and the performance of the high order nonlinear filter schemes with pre- 
processing and post-processing steps in conjunction with the use of a high order non- 
dissipative spatial base scheme [51, 54] are briefly illustrated for two test cases. 

2.1. High Order Nonlinear Filter Schemes [34, 49, 51, 54] 

Before the application of a high order non-dissipative spatial base scheme, the pre- 
processing step to improve stability had split inviscid flux derivatives of the governing 
equation(s) in the following three ways, depending on the flow types and the desire for 
rigorous mathematical analysis or physical argument. 

• Entropy splitting of Olsson & Oliger [25] and Yee et al. [45, 46]: The resulting 
form is non-conservative and the derivation is based on entropy norm stability with 
numerical boundary closure for the initial value boundary problem. 
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• The system form of the Ducros et al. splitting [6]: This is a conservative splitting and 
the derivation is based on physical arguments. 

• Tadmor entropy conservation formulation for systems (Sjogreen & Yee [35]): The 
derivation is based on mathematical analysis. It is a generalization of Tadmor’s en- 
tropy formulation to systems and has not been fully tested on complex flows. 

After the application of a non-dissipative high order spatial base scheme on the split 
form of the governing equation(s), to further improve nonlinear stability from the non- 
dissipative spatial base scheme, the post-processing step of Yee & Sjogreen [49, 51], 
Sjogreen & Yee [34] nonlinearly filtered the solution by a dissipative portion of a high order 
shock-capturing scheme with a local flow sensor. These flow sensors provide locations and 
amounts of built-in shock-capturing dissipation that can be further reduced or eliminated. 
For all the computations shown, the Ducros et al. splitting is employed since a conserva- 
tive splitting is more appropriate if one does not know if the subject flow is shock-free or 
turbulence with shocks. Some attributes of the high order filter approach are: 

• Spatial Base Scheme: High order and conservative with high order freestream preser- 
vation metric evaluation for curvilinear grids, (no flux limiter or Riemann solver) 

• Physical Viscosity: Automatically taken into consideration by the base scheme. The 
same order of central differencing for the viscous derivative as the convective flux 
derivatives are used. 

• Efficiency: One Riemann solve per dimension per time step, independent of time dis- 
cretizations (less CPU time and fewer grid points than their standard shock-capturing 
scheme counterparts) 

• Accuracy: Containment of numerical dissipation via local wavelet flow sensor 

• Well-balanced scheme: These nonlinear filter schemes are well-balanced schemes for 
certain chemical reacting flows and problem containing geometric source terms [39] 

• Stiff Combustion with Discontinuities: For some stiff reacting flow test cases, the 
high order filter scheme is able to obtain the correct propagation speed of discontinu- 
ities whereas the standard high order WENO scheme cannot [19, 56]. 

• Parallel Algorithm: Suitable for most current supercomputer architectures 

2.2. Sample test Cases Illustrating the Efficiency and Accuracy of High Order Filter 
Schemes 

These filter schemes are efficient, and the total computational cost for a given error tol- 
erance is lower than for standard shock-capturing schemes of the same order. This is of 
importance, for example, in DNS and in flow control optimization to improve aerodynamic 
properties, where the flow simulation must be carried out many times during the optimiza- 
tion loop. The efficiency and accuracy of the schemes for a wide variety of flow problems 
can be found in aforementioned cited references. Here two test cases are illustrated. 
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Figure 1: CPU comparison of four shock-capturing schemes. 


2.2.1. 2D Shock/Vorticity Interaction 

Figure 1 shows a comparison of a second-order TVD, seventh-order WENO (WEN07), 
hybrid scheme (switch between eighth-order spatial central scheme and WEN07 using 
wavelet flow sensor as the switch indicator) and the filter scheme WEN07fi (an eighth- 
order spatial central base scheme and the dissipative portion of WEN07 and the same 
wavelet flow sensor to guide where the WEN07 dissipation should be applied at the post- 
possessing nonlinear filter step). A second-order Runge-Kutta method was used for the 
TVD scheme and the classical fourth-order Runge Kutta method was used for the rest of 
the spatial scheme. For this particular simple 2D shock-vorticity interaction test case with a 
simple weak planar shock without structure, WEN07, hybrid, and WEN07fi give the same 
accuracy. However, there is large gain in CPU time by the filter scheme for this turbulence- 
free test case. For turbulence with shocks, there is a more beneficial gain both in accuracy 
and CPU time of the filter schemes over the their standard WENO counterparts. 

2.2.2. ID Shock/Turbulence Interaction Problem 

This 1-D compressible inviscid ideal gas problem is one of the most computed 
test cases in the literature to assess the capability of a shock-capturing scheme in the 
presence of shock/turbulence interactions. The flow consists of a shock at Mach 3 
propagating into a sinusoidal density field with initial data given by (p L , u L , p L ) = 
(3.857143, 2.629369, 10.33333) to the left of a shock located at x = —4, and 
( pr , ur, pr ) = (1 + 0.2sin(5x), 0, 1) to the right of the shock, where p is the den- 
sity, u is the velocity and p is the pressure. The computational domain is [—5, 5] and the 
computation stops at time equal to 1.8. Figure 2 shows the comparison among WEN03, 
WEN05 and WEN07, and their corresponding filter schemes WEN03fi, WEN05fi and 
WEN07fi using a very coarse uniform grid of 200 points with the reference solution. The 
reference solution is obtained with WEN05 using 16000 grid points. WEN05fi required 
at the most 50% of the CPU time of WEN05 if third or fourth-order Runge-Kutta time 
discretization were used. In order for WEN05 to obtain a similar accuracy as WEN05fi, at 
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Figure 2: 


least two times the number of grid points is needed. Moreover, the accuracy of WEN05fi 
is similar to WEN09 (computation not shown). 

2.3. Objective and Outline 

The objective of this paper is to use the same TML problem setup and convective 
Mach cases as in [13] to compare the performance of WEN07fi with standard WEN05 
and WEN07 for convective Mach M c = 0 . 1 , 0 . 3 , 0 . 8 , 1 . 0 , 1 . 5 . For WEN07H, no tuning of 
scheme parameters is needed for all theM c cases. For detailed physics, see [13]. 

The outline of this paper is as follows: The high order nonlinear filter methods are sum- 
marized in Section 2. Recent improvement of the scheme will also be briefly discussed. 
The problem setup for the temporally resolving mixing layer (TML) is given in Section 
3. Numerical results comparing the performance among WEN05, WEN07 and the as- 
sociated WEN07fi are then presented in Section 4. The comparison is focused only on 
the momentum thickness and compressibility factor as a function of the studied convective 
Mach numbers. Even WEN05 is two order lower than WEN07 WEN07fi, it is chosen for 
the comparison as a benchmark study. This is due to the fact that WEN05 has been widely 
used for many applications. 

3. Numerical Methods 

This section summarizes the numerical methods to be used for the turbulent TML study. 
The numerical methods solve the split form of the inviscid flux derivatives according to the 
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pre-processing step. The discussion is broken up into two subsections. 
3.1 . Original High Order Filter Method 


For turbulence with shocks, instead of solely relying on very high order high-resolution 
shock-capturing methods for accuracy, the filter schemes [44, 45, 34, 48, 49] take advantage 
of the effectiveness of the nonlinear dissipation contained in good shock-capturing schemes 
as stabilizing mechanisms at locations where needed. Such a filter method consists of two 
steps: a full time step using a spatially high-order non-dissipative base scheme, followed 
by a post-processing filter step. The post-processing filter step consists of the products 
of wavelet-based flow sensors and nonlinear numerical dissipations. The flow sensor is 
used in an adaptive procedure to analyze the computed flow data and indicate the location 
and type of built-in numerical dissipation that can be eliminated or further reduced. The 
nonlinear dissipative portion of a high-resolution shock-capturing scheme can be any TVD, 
MUSCL, ENO, or WENO scheme. By design, the flow sensors, spatial base schemes 
and nonlinear dissipation models are standalone modules. Therefore, a whole class of low 
dissipative high order schemes can be derived with ease. Unlike standard shock-capturing 
and/or hybrid shock-capturing methods, the nonlinear filter method requires one Riemann 
solve per dimension, independent of time discretizations. The nonlinear filter method is 
more efficient than its shock-capturing method counterparts employing the same order of 
the respective methods. 

Recently, these filter schemes were proven to be well-balanced schemes [39] in the 
sense that these schemes preserve exactly certain steady state solutions of the chemical 
nonequilibrium governing equation. With this added property these filter schemes can 
better minimize spurious numerics in reacting flows containing mixed steady shocks and 
unsteady turbulence with shocklet components than standard non-well-balanced shock- 
capturing schemes. In addition, for some stiff reacting flow test cases, the high order filter 
scheme is able to obtain the correct propagation speed of discontinuities whereas the stan- 
dard high order WENO scheme cannot [19, 56]. 

For simplicity of the presentation the discussion for the base scheme and post- 
processing step of the filter scheme is restricted to the inviscid part of the Navier-Stokes 
equations. For viscous gas dynamics the same order of spatial centered base scheme for the 
convection terms and the viscous terms are employed. For all of the LES computations the 
classical fourth-order Runge-Kutta time discretization is employed. 


Consider the 3-D compressible Euler equations in Cartesian geometry, 


U t + V • F = 0; U 



( pu \ 

pnn T + p . 
\u(e +p) / 


( 1 ) 


Here the velocity vector u = (u, v , w) T , the momentum vector m = (pu, pv, pw), p is the 
density, and e is the total energy. 
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In a Cartesian grid denote the grid indices for the three spatial directions as (j, k, l). 
The spatial base scheme to approximate the x inviscid flux derivatives F(U) X (with the grid 
indices k and l for the y- and ^-directions suppressed) is written as, e.g.. 


dF 

dx 


Dq&Fj, 


( 2 ) 


where D 0 8 is the standard eighth-order accurate centered difference operator. See [35] for 
the split form of 2. 

After the completion of a full Runge-Kutta time step of the base scheme step, the second 
step is to adaptively apply a post-processing nonlinear filter. The nonlinear filter can be 
obtained e.g., in the ^-direction by taking the full seventh-order WENO scheme (WEN07) 
[16] for the inviscid flux derivative in the £ -direction and subtracting D 08 Fj. The final 
update of the solution is (with the numerical fluxes in the y- and ^-directions suppressed as 
well as their corresponding y- and ^-directions indices on the x inviscid flux suppressed) 


At 


= izW+1,2 


H* 


( 3 ) 


Ax^ 3+l/z i- 1 / 2 *' 

The nonlinear filter numerical fluxes usually involve the use of field-by-field approxi- 
mate Riemann solvers. If the Roe type of approximate Riemann solver [30] is employed, 
for example, the a -filtcr numerical flux vector Hj + 1 / 2 evaluated at the U* solution from the 
base scheme step is 


Hj+ 1/2 = Rj+l/2Hj+i/2, 

where Rj+ 1/2 is the matrix of right eigenvectors of the Jacobian of the inviscid flux vector in 
terms of the Roe’s average states. Denote the elements of the vector H J+l / 2 by //■ H / 2 , l = 
1, 2, ..., 5. The nonlinear portion of the filter ^+ 1/2 l ias the form 

h j+ 1/2 = 2 a; 5+V2^+l/2- ( 4 ) 

Here u l - + 1 , 2 is the wavelet flow sensor to activate the nonlinear numerical dissipation 
\4 >l j+ 1/2 an4 the original formulation for n is a positive parameter that is less than or equal 
to one. Some tuning of the parameter k is needed for different flow types. A local n to be 
discussed next, depending on the local Mach number for low speed flows and depending 
on local shock strength for high speed flows, would minimize the tuning of parameters. A 
local flow sensor was discussed by Lo et al. [24] by taking advantge of the Ducros et al. 
shock flow sensor [7] to obtain a local artificial compression method (ACM) sensor for the 
original Yee et al. filter scheme [44] . 

The dissipative portion of the nonlinear filter \4> l j + \/2 = fij+ 1/2 — ^j+ 1 / 2 ' s the dissipative 
portion of, e.g., WEN07 for the local /th-characteristic wave. Here fj' l+l /2 an <3 ^j+1/2 are 
numerical fluxes of WEN07 and the eighth-order central scheme for the /th characteristic, 
respectively. Hereafter, we denote this filter scheme as WEN07fi. 
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Mach sensors, e=0 



Figure 3: Mach number sensors. f(M) (blue) function by Li and Gu, fi(M) (red) modified /(M), and 
/ 2 (M) (black) (includes low supersonic Mach numbers). 


A summary of the three basic steps to obtain o/ +]/ , 2 can be found in [34, 49]. For ex- 
ample, the flow sensor u;*. +1 , 2 to activate the shock-capturing dissipation using the cut off 
procedure is a vector (if applied dimension-by-dimension) consisting of “l’s” and “0’s”. 
For all of the computations, a three-level second-order Harten multiresolution wavelet de- 
composition of the computed density and pressure is used as the flow sensor [34]. 

3.2. Improved High Order Filter Method 

Previous numerical experiments on a wide range of flow conditions [44, 45, 34, 48, 49] 
indicated that the original filter scheme improves the overall accuracy of the computation 
compared with standard shock-capturing schemes of the same order. Studies found that 
the improved accuracy is more pronounced if the parameter k in (4) is tuned according 
to the flow type locally. For hypersonic flows with strong shocks, k is set to 1. For high 
subsonic and supersonic flows with strong shocks, k is in the range of (0.3, 0.9). For low 
speed turbulent flows without shocks or long time integration of smooth flows, k can be 
one to two orders of magnitude smaller than 1. In other words, k should be flow location, 
shock strength and local flow type dependent. The improved k proposed in [51] consists of 
a simple global k for smooth flows and a local k for problems with shocks and turbulence. 

3.2.1. An efficient global nfor low Mach number and smooth flows: 

The flow speed indicator formula of Li & Gu to overcome the shortcomings of “low 
speed Roe scheme” [23] was modified to obtain an improved global k denoted by 7c for (4) 
to minimize the tuning of the original n for low Mach number flows. 7c has the form: 
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with 


« = /i(M)k, 


(5) 


/i(M) 


min 


/M 2 y/4 + (1 - M 2 ) 2 
1 + M 2 


(6) 


Here M is the maximum Mach number of the entire computational domain at each stage 
of the time evolution. fi(M) has the same form as [23] except there is an extra factor “4p’ 
added to the first argument on the right-hand-side of the original form f(M) in equation 
(18) of [23]. The added factor provides a similar value of the tuning k observed from nu- 
merical experimentation reported in aforementioned cited references. With the flow speed 
indicator f\ (M) in front of k, the same n used for the supersonic shock problem can be used 
without any tuning for the very low speed turbulent flow cases. Another minor modification 
of the above is, 


fi(M) = max 


min 


/M 2 yfl + (1 — M 2 ) 2 
\~2 1 + M 2 


5 


where e is a small threshold value to avoid completely switching off the dissipation. A 
function which retains the majority of f\ (M) but includes larger Mach number for not very 
strong shocks is 

/ 2 (M) = (Q(M,2) + Q(M,3.5))/2 


or 

h W) = ma x((Q(M, 2) + Q(M, 3.5))/2, e), 


where 


The polynomial 


Q(M, a ) 


P(M/a) M < a 

1 otherwise 


P(x) = x 4 (35 - 84a; + 70a; 2 - 20a; 3 ) 


is monotonically increasing from P( 0) = 0 to P( 1) = 1 and has the property that P'(x) 
has three continuous derivatives at x = 0 and at x = 1. 


Below supersonic speeds, a simple and efficient global k can be obtained according to 
the maximum Mach number of the entire flow field and the value is determined by f\ (M) 
or / 2 (M) for non-zero a/ , |/2 . It is noted that if the original f(M) were used instead of 
fi(M) or / 2 (M) in Eq.(5), the amount of nonlinear filter dissipation could be too large for 
very low speed turbulent flows (for the same fixed k). See Fig. 3 for details. 
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3.2.2. Local flow sensor for a wide spectrum of flow speed and shock strength 

At each time step and grid point, the aforementioned global 7c is not sufficient to reduce 
the amount of numerical dissipation where needed for flows that contains a variety of flow 
features. A more appropriate approach is to obtain a“local k ” that is determined according 
to the above at each grid point. If known, a dominating shock jump variable should be used 
for shock detections. In other words, the filter numerical flux indicated in Eq.(4) is replaced 
by: 


h j+ 1/2 


9 K + 1/2^ + 1/2^+1/2]- 


(7) 


In the case of unknown physics and without experimental data or theory for compari- 
son, k l j + i/ 2 has to depend on the local Mach number in low speed or smooth flow regions, 
depend on local shock strength in shock regions and depend on turbulent fluctuations in 
vortical regions in order to minimize the tuning of parameters. According to the flow type 
locally, for each non-zero wavelet indicator uJ l j+1 / 2 , ^+ 1/2 should provide the aforemen- 
tioned amount (between (0, 1)) to be filtered by the shock-capturing dissipation <t> l j+1 / 2 - F° r 
problems containing turbulence and strong shocks, the shock strength should come into 
play. One measure of the shock strength can be based on the numerical Schlieren formula 
[12] for the chosen variables that exhibit the strongest shock strength. In the vicinity of 
turbulent fluctuation locations, n l - +l , 2 will be kept to the same order as in the nearly incom- 
pressible case, except in the vicinity of high shear and shocklets. 


Due to the fact that k works well for local Mach number below 0.4, k only needs to be 
modified in regions that are above 0.4. In other words, the final value of n l - +l , 2 is determined 
by the previous local K, if the local Mach number is below 0.4. If the local Mach number 
is above 0.4, at discontinuities detected by the non-zero wavelet indicator uj l j +1 , 2 , k 7 +1/2 
is determined by the shock strength (normalized between (0, 1)) based on the Schlieren 
formula near discontinuities. Again, if known, dominating shock jump variables should be 
used for shock detections. At locations with turbulence, determined by the turbulent sensor 
(e.g., u l ]+l j 2 obtained from employing wavelets with higher order vanishing moments), 
K l j + 1 / 2 is kept to the same order as in the nearly incompressible case, except in the vicinity 
of high shear and shocklet locations, where a slightly larger n l j+ i/ 2 would be used. Methods 
in detecting turbulent flow can be (a) Wavelets with higher order vanishing moments, and 
(b) Wavelet based Coherent Vortex Extraction (CVE) of Farge et al. [8] (Split the flow into 
two parts: Active coherent vortices and incoherent background flow). 

Results by the local k 1 j +1 , 2 that take the local flow speed and shock strength into consid- 
eration will be reported in [53], an expanded version of [51]. Preliminary study with more 
complex shock turbulence problems and the applicability of even wider flow types indicates 
the necessity of the local K l j+1 / 2 . 

In this paper, all the computations use the global K, the Ducros et al. splitting of the 
inviscid flux derivatives and WEN07fi using the global Ti in conjunction with the classical 
fourth-order Runge-Kutta temporal discretization, k = 0.7 for all test cases. 
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4. Description of the Physical Problem and Computational Setup 

LES using the subgrid model of [27] of a TML between two streams with equal and 
opposite velocities was considered in [13]. In [13] the three main characteristics of com- 
pressible TML (the self similarity property, compressibility effects and the presence of 
large-scale structure with shocklets for high M c ) were considered for the LES study. The 
role of compressibility in turbulent mixing layers remains an important issue in aeronautics. 
Lor example, in the design of high-speed (supersonic or hypersonic) propulsion devices, the 
stabilizing effect of compressibility may reduce the efficiency of engines in mixing the fuel 
and the oxidizer. One of the objectives of the study in [13] is to use WEN07fi to in- 
vestigate the compressibility effects in highly sheared turbulent flows subjected to strong 
shocks. Here WEN07fi refers to the pre-processing step (Ducros et al. splitting of the 
inviscid flux derivative) in conjunction with the eighth-order central spatial base scheme 
with the dissipative portion of WEN07 and the global flow sensor discussed in Section 2 
as the post-processing nonlinear filter step. The objective of the current investigation is to 
compare the performance among WEN07fi, WEN05 and WEN07 using the same problem 
setup, computational domain and grid size as in [13]. 

4.1. TML Problem Setup in [13 ] 

The configuration of the TML is shown in Lig.4. Live test cases are carried out with 
different convective Mach numbers (M c = 0.1, 0.3, 0.8, 1.0 and 1.5) ranging from the in- 
compressible case M c = 0.1 up to the supersonic case M c = 1.5. The M c = 1.5 case 
correspond to a highly compressible mixing layer, whereas the first two cases M c = 0.1 
and M c = 0.3 can be considered as quasi-incompressible and are used to compare with 
the experimental results of an incompressible shear layer. All of the simulations described 
below are performed at an initial Reynolds number, Re wo , based on the mean velocity dif- 
ference A U = 2Ui = — 2f/ 2 , the average viscosity of the free streams and the vorticity 
thickness <5 Uo of 800 with S Wo = 4 6g 0 , where A. = AU / (du/dy) max is the vorticity thick- 
ness of the shear layer, and 5g is the momentum thickness (see [13] for details). The values 
of Re w can be as large as 3 x 10 5 at the end of the simulation, which is one order of magni- 
tude higher than similar DNS and LES computations reported in the literature [28, 26, 9]. 
The mean flow is initialized with a tangent hyperbolic profile for the streamwise velocity, 
u(y ) = \AU tanh [y/(2Sg 0 )\, while the two other velocity components are set to zero. In 
addition to these mean values, three-dimensional turbulent fluctuations (u',v',w') are im- 
posed, while initial pressure and density are set constant. Since the simulation is temporal, 
the initial velocity perturbations are computed using a digital filter technique [18]. This 
procedure utilizes the prescribed Reynolds stress tensor and length scales of the problem 
concerned to generate the corresponding fluctuating velocity field, taking into account the 
nature of the autocorrelation function for the prevailing turbulence. See [13] for details. The 
digital filter algorithm is given in Appendix B of [13]. The length scales are chosen as S Wo 
in each direction. The Reynolds stresses have a Gaussian shape in y with amplitudes chosen 
to be similar to the experimental peak intensities observed in incompressible mixing layer 
[ 2 ], 
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Figure 4: Schematic configuration of temporal mixing layer. 

Periodic boundary conditions are enforced in the streamwise (;/;) and spanwise (z) di- 
rections, while non-reflecting conditions are applied at both top and bottom boundaries ( y 
direction). The use of a periodic boundary condition in the x direction corresponds to the 
temporal formulation of mixing layer evolution, which evolves only in time as it spreads in 

y- 

4.2. Mesh Requirements 

Similar to [9], a computational domain of lengths L x x L y x L z = 1200 5g 0 x 370 Sg 0 x 
270 Sg 0 is used with the corresponding mesh points N x x N y x N z = 512 x 211 x 131. 
The same grid system uniformly spaced in the x and x directions and stretched in the y 
direction is employed for all considered cases. The stretching function for the y-dircction 
is based on ■ where L y is the box size in the ^/-direction and the stretching factor 

by = 3.4. The mapped coordinate // is equally spaced and runs from —1 to 1. The grid used 
in this study contains an order of magnitude fewer cells than that of the DNS of Pantano and 
Sarkar [28] compared to the domain length. To ensure that the computational domain in the 
x- and z-dircctions is sufficiently wide, the two-point correlation functions were analyzed 
in [13]. 

5. Numerical Results 

LES using the subgrid model of [27] and WEN07fi was performed in [13] on the TML 
problems at different convective Mach numbers M c (0.1, 0.3, 0.8, 1.0, 1.5). Studies in [13] 
show the level of good agreement obtained between LES and DNS for convective Mach 
number 0.3. In higher convective Mach numbers (up to 1.5), LES results are in good agree- 
ment with the experiments. Also, the principal compressibility effects such as the reduction 
of the spreading rate and the turbulence intensities are well predicted. 
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Figure 5: 2-D cut at midplane of instantaneous dilatation flow-field at r = 1000 for three different convective 
Mach numbers (M c = 1.5 (top), 1.0 (middle), and 0.8 (bottom)). 
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The flow features of this TML are determined by M c . LES computations are carried 
out up to dimensionless time r = tAU / Sq 0 — 3000 for the higher convective Mach number 
cases and r ~ 1200 for the quasi-incompressible cases. Shocklets developed for M c > 0.8. 
The shocklets are stronger and more complicated as M c increases. Figure 4 shows the 2-D 
cut at midplane of instantaneous dilatation flow field at r = 1000 for three different convec- 
tive Mach numbers (M c = 1.5, 1.0, and 0.8). Figure 4 indicates the different flow pattern 
as a function of M c . Figure 5 shows the 2-D cut at mid plane and 3-D of instantaneous nu- 
merical Schlieren of vorticity at r = 2000 and M c = 1.5 by WEN07A. Figure 6 shows 3-D 
iso-surfaces of instantaneous vorticity field from different viewing angles for M c = 1.5. 

For each M c , after a transient phase, the momentum thickness approaches a separate 
linear growth regime. FES results by WEN07fi for different convective Mach numbers 
agree well with the analytically predicted slopes whereas WEN05 and WEN07 do not 
agree well with the predicted slope for most of the M c cases. 

Figures 7-11 show the momentum thickness comparison among the three schemes for 
each of the M c cases. It is interesting to note that the global k used by WEN07A only shows 
a slight improvement over WEN05 and WEN07 for M c = 0.8 and 1.0. For the compress- 
ibility factor computations to be shown next, WEN07A compares well with experimental 
data for all studied M c , whereas this is not the case for WEN05 and WEN07. 

Figure 12 displays the compressibility factor as a function of M c by three high order 
schemes comparing with published work and experiments. This figure shows the superior 
performance of WEN07fi compared with WEN05 and WEN07. The FES results using 
WEN07fi agree well with experimental results of Barone et al. (2006) [1], and published 
direct numerical simulations (DNS) work of Rogers & Moser (1994) [31] and Pantano & 
Sarkar (2002) [28]. In all M c cases, no tuning of WEN07A scheme parameters was needed. 
For all the M c cases considered, solutions by WEN05 and WEN07 compared poorly with 
experimental data and DNS computations. 

The unsteady time evolution of turbulence with shocklets for M c = 0.8, 1.0 and 1.5 
among WEN05, WEN07 and WEN07fi are very different in terms of the location and 
strength of the shocklets and turbulent fluctuation pattern. Figure 13 shows the 2-D cut 
at midplane of the numerical Schlieren of vorticity at r = 2000 and M c = 1.5 computed 
by WEN07A and WEN05. Figure 14 shows the 2-D cut at midplane of instantaneous 
dilatation flow-field at r = 2000 and M c = 1.5 computed by WEN07A and WEN05. 

6. Concluding Remarks 

The present work serves as a validation and performance study of the improved filter 
schemes of [51] on a representative complex compressible turbulent flow consisting of a 
wide range of flow speeds. All the computations use the global 7 1, the Ducros et al. splitting 
of the inviscid flux derivatives and WEN07A with k = 0.7 described in Section 2.2.1. In all 
M c cases, no tuning of WEN07A scheme parameters were needed. Studies indicated that 
WEN07A compared well with experimental data and published DNS work. For all the M c 
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wave angle 



Figure 6: Instantaneous numerical Schlieren pictures at r = 2000 and M c = 1.5 by WEN07fi indicating 
shocklets formation. 
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Figure 7: Temporal mixing: Instantaneous numerical Schlieren of vorticity by WEN07fi at r = 2000, M c = 
1.5, top and side views. 
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Figure 8: Temporal mixing: Momentum thickness comparison for M c = 0.1 



Figure 9: Temporal mixing: Momentum thickness comparison for M c = 0.3 
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Figure 10: Temporal mixing: Momentum thickness comparison for M c = 0.8 



Figure 11: Temporal mixing: Momentum thickness comparison for M c = 1.0 
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Figure 12: Temporal mixing: Momentum thickness comparison for M c = 1.5 


cases considered, solutions by WEN05 and WEN07 compared poorly with experimental 
data and DNS computations. 

The same high order filter scheme is being used for the simulation of two much higher 
M c cases of M c = 2, 3. The computational box size, especially in the ^/-direction has to 
be doubled or more. A finer grid is also needed in order to obtain an accurate and stable 
solution. These computations are many times more CPU-intensive than the lower M c cases. 
Results will be reported in a forthcoming paper. 
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